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ABSTRACT 

We present a pulsar-atmosphere (PA) model for modulated thermal X-ray emission from 
cooling magnetized neutron stars. The model synthesizes the spectral properties of detailed 
stellar atmosphere calculations with the non-uniform surface properties anticipated for isolated, 
aging radio pulsars, and general relativistic effects on photon trajectories. We analyze the archival 
Chandra observations of the middle-aged radio pulsar PSR 0656-1-14 with the PA model and find 
it is an excellent representation of the phase-averaged X-ray spectrum for a broad range of polar 
effective temperature Tp and column density Nh — 10^° cm~^. The spectral fits favor a sub- 
solar neutron star mass M < 1.0 Mq, a large radius R — 15—16 km for a source distance 
d ^ 190 pc, while retaining consistency with theoretical neutron star equations-of-state. The 
modulated spectrum constrains the angular displacement of the pulsar spin axis to 30 ± 2° with 
respect to the line of sight. 



1. Introduction 

The discovery of modulated thermal X-rays in ROSAT observations of a handful of neutron stars (NSs) 
provided early evidence that the surface properties of these sources are inhomogeneous. Most NSs are 
rotation powered, and emit predominantly non-thermal X-rays generated exterior to the star in the pulsar 
magnetosphere (Becker & Triimper 1997; Possenti et al. 2002). Thermal X-rays from a few exceptional 
sources originate at the stellar surface and can be distinguished from magnetospheric contributions owing 
to their intermediate age: these are sufficiently aged to have experienced a decline in synchrotron activity, 
yet not so old that they are altogether too faint to be seen. The "Three Musketeers," PSR 0656-1-14, PSR 
1055-52 and PSR 0630+18 (Geminga) are relatively bright and weakly absorbed middle-aged radio pulsars 
which are among the earliest known thermal emitters (see Ogelman 1995 for a summary of ROSAT data). 
These sources remain excellent candidates for evaluating NS cooling models and surface characteristics, and 
are objects of continued interest in Chandra and XMM-Newton observations. Becker & Triimper (1997) 
tabulated the short catalog of confirmed thermal NSs from the combined ROSAT and ASCA data sets, 
which include the relatively young Vela pulsar. These X-ray spectra are found to be poorly described 
by a single blackbody and require in addition either a power-law or second blackbody component. The 
secondary "hard" thermal component is invoked to bridge the energy range between the "soft" thermal 
and power-law contributions and broadens the spectral distribution of the soft X-rays. This spectral model 
may be interpreted as thermal radiation from the entire stellar surface (soft) seen in relief against small 
hot-spots or otherwise externally heated polar regions (hard) , plus the contribution from the magnetosphere 
at sufficiently high energies. In this picture, a careful measurement of the soft thermal component may allow 
for direct measurement of intrinsic stellar properties - radius R, characteristic temperature T, magnetic field 
and surface composition. This is quite challenging in practice because a reasonable interpretation of the 
spectrum requires a self-consistent model of the pulsar's synchrotron continuum, stellar viewing angles and 
distance, in addition to the NS surface properties and the associated thermal radiation. 
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The principal difficulty posed by the soft X-ray spectra of these pulsars is in reconciling their modulated 
thermal emission and estimated distance with a model for the surface radiation. Owing in part to modest 
instrumental sensitivity to the peak energies of their relatively cool surface radiation (0.1 — 1.0 keV), measures 
of temperature, the intervening absorbing column and overall normalization (i.e., radius and distance) are 
mutually affected. These uncertainties can be compounded by discrepancies as large as a factor of ^ 4 from 
different techniques for estimating the distance to some sources. The average spectrum alone is not sufficient 
to deduce the temperature at the stellar surface, where local anisotropics likely play a significant role in the 
fraction and profile of the pulsed radiation; the energy dependent pulse profiles and modulation index are 
functions of the apparent viewing geometry and radiation beaming patterns, but are also dependent on the 
assumptions of the model for the relevant pulsar mechanism. In the absence of perfect energy resolution 
of the detector, the strength of the modulation that the observer infers also depends on the amount of 
interstellar absorption (Page 1995; Perna et al. 2000). 

Models with two distinct thermal components provide acceptable fits to the mean X-ray flux of these 
sources, while more complicated surface temperature distributions also fare well. Page (1995) considered a 
smoothly varying blackbody distribution on the stellar surface and noted that a non-uniform temperature 
distribution implies increasing pulsed fraction with energy. The integration of such a dipolar temperature 
blackbody distribution is well represented by that of a single temperature within the limits of the PSPC 
data sets. X-ray pulsations are suppressed by gravitational lensing of the stellar surface for a dipolar field 
regardless of the details of radiation beaming. More complicated field geometries can be invoked to explain 
observed X-ray pulsation from the radio pulsars, but require unacceptable tuning of surface parameters (Page 
& Sarmiento 1996). The blackbody temperature configurations generate a pulsed fraction which increases 
with energy, and seems incapable of explaining either the non-monotonic growth in pulsed flux seen for the 
middle-aged radio pulsars observed with ROSAT and Chandra, or the phase shift in pulse synchronization 
seen at 0.3 — 0.5 keV in ROSAT observations. These may be partially resolved by beaming of radiation 
owing to intrinsic anisotropics to explain both phase-resolved measures. 

Radiation emergent from a stellar atmosphere is inherently anisotropic regardless of the local tempera- 
ture and magnetic field and, when combined with a non-uniform temperature proffie, is capable of generating 
considerable pulsar modulation. The earliest effort to describe X-ray pulsations from NSs with magnetized 
plasma atmospheres was due to Pavlov et al. (1994). These models characterized the pulsar emission with 

small polar caps having uniform normal magnetic field while the remainder of the stellar surface was assumed 
to be too cool to contribute to the total flux. Zavlin et al. (1995) refined this approach and, in particular, 
generalized their calculation to include an inhomogeneous surface temperature distribution and dipolar field 
geometry. A significant result of this effort was the discovery of a "crossover" energy at which the pulsed 
X-ray spectrum could experience a phase-shift; this transition depends critically upon the details of the ion- 
ization equilibrium of cool NS atmospheres, and on the assumed surface temperature profiles and magnetic 
field geometry. The dipole model was found to be generally consistent with the ROSAT observations of the 
middle-aged radio pulsars. Meyer et al. (1994) made a detailed fit to the ROSAT observations of Geminga 
using partially ionized magnetic hydrogen atmosphere models; these aiithors considered a simple model for 
Geminga in which the temperature and magnetic field of the star were assumed homogeneous, effectively a 
magnetic monopole. The harder spectrum of the plasma atmosphere was not alone able to fully describe 
the ROSAT PSPC data and required the secondary polar cap component, albeit with much smaller effective 
emitting area than that deduced from dual blackbody fits. 

Perna et al. (2001) derived pulsed X-ray spectra (0.1 — 1.0 keV) of rotating NSs from the surface 
temperature profiles of Heyl & Hernquist (1998a) by integrating the spectrum of a toy hydrogen atmosphere 
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model (Heyl & Hernquist 1998b) at each point on the stellar surface. As found in the blackbody analysis, two 
thermal components were required to reproduce the averaged ROSAT spectra of the radio pulsars; the phase 
resolved properties of these spectra (pulsed fraction and pulse phase-lag with energy) were also reproduced 
by assuming the cap had different beaming properties than the remainder of the surface. More specifically, 
Perna et al. (2001) foimd that the behavior of the pulsed fraction, as measured by ROSAT PSPC, could 
be reproduced with the combination of a "pencil beam" component generated over the entire surface of the 
star, and a hotter "fan" beamed component, produced at the polar caps. 

1.1. PSR 0656+14 

X-ray emission from PSR 0656+14 was first detected with EXOSAT (Cordova et al. 1989); the suggested 
X-ray pulsations were later confirmed when evidence of X-rays pulsed at the radio period P ~ 385 ms was 
found in the ROSAT PSPC observations described by Finley et al. (1992); this spectrum was inferred to be 
consistent with a blackbody T^i, ~ 9 x 10^ K having normalization consistent with R = IQ km assuming a 
source distance d = 500 pc. An optical counterpart was found in ESQ observations {my 25) as reported by 
Caraveo et al. (1994), and optical pulsations displaced from the radio phase by 0.2 were recovered from later 
observations (Shearer et al. 1997). At X-ray energies, the correlation between spectral hardness and rotation 
phase suggests the pulsation originates in a non-uniform surface temperature distribution. Anderson et al. 
(1993) compare these PSPC data to contemporary ROSAT HRI observations (0.15 — 2.0 keV) and conclude 
from the modulated X-rays that the pulsed emission is probably no harder than the unpulsed spectrum. 
Greiveldinger et al. (1996) considered joint fits to the ROSAT and ASCA data sets found that the spectrum 
was best fit by two thermal components ~ 8 x 10^ K, Th — 1.5 x 10^ K for Nh.20 = 1-7 (expressed in units 
of 10^" cm^^); the fit was significantly improved by inclusion of a power-law to accommodate the hard X-ray 
tail seen with ASCA. Edelstein et al. (2000) find that simultaneous fits to EUVE and ROSAT observations 
require both cooler temperatures for the soft thermal component Tbt — 3.6 x 10^ K and larger absorbing 
columns Nh,20 = 2.21; continuing the fit to the thermal optical measurements (Pavlov et al. 1996, 1997) 
allow for larger Tjj = (4.7 — 7.4) x 10^ K and smaller Nh,20 = (1-1 — 1-8), consistent with Greiveldinger et al. 
(1996). In either case, the hard thermal component is ~ (1.1 - 1.2) X 10^ K. Chandra LETG observations 
(Marshall & Schulz 2002) find double blackbody fit parameters consistent with those of Greiveldinger et al. 
(1996). Optical photometry performed by Kurt et al. (1998) suggests the non-thermal optical flux has a 
spectral index roughly consistent with that of the hard X-ray emission; subsequent IR, optical and near UV 
measurements are consistent with a harder spectrum and evolution of the spectral index (Koptsevich et al. 
2001). Pavlov et al. (2002) also made an ACIS-S observation in CC mode, and examined the joint fit with 
the earlier grating data (0.2 — 6.0) keV. The resultant three component (BB+BB+PL) model is likewise 
consistent with that of Greiveldinger et al.; these authors find a fit to dual magnetic atmosphere models 
questionable because of the large derived stellar radius, assuming the source distance d ~ 500 pc. HST 
observations of PSR 0656+14 provide strong evidence for identification of the optical counterpart (Mignani 
et al. 2000). The proper motion of this counterpart has been measured but, owing to the uncertainty in 
the pulsar velocity, does not constrain the source distance to better than a factor of ~ 3. The lower bound 
likewise agrees with the photometric estimate for this pulsar (150 — 200 pc) by Golden & Shearer (1999). 

The energy dependence of the pulsed fraction of thermal X-rays has been measured from Chandra 
ACIS-S observations in the continuous clocking mode (Pavlov et al. 2002). The pulsed fraction does not 
grow continuously with energy, but rather exhibits lower pulsation from 0.5-1.0 keV when compared to 
either lower or higher energy intervals; this property is a significant constraint on comprehensive models for 
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the X-ray emission. The earlier ROSAT observations show an increase in the pulsed fraction with energy, 
and is seemingly inconsistent with any decrease in pulsation on this same interval. Both datasets reveal an 
apparent phase shift in the pulsed thermal X-rays although the magnitude of the phase displacement differs 
and energies affected differ between the two observations 0.8 keV in ROSAT, > 1.0 keV in Chandra). 
Possenti ct al. (1996) evaluated the pulsed fraction predicted by the dual temperature blackbody model and 
concluded from the PSPC data that the measured ~ 9% on (0.1 — 2.4) keV could be reconciled with the 
viewing geometry x = 30 estimated by Rankin (1993), but only for gravitational redshift z < 0.15. 

In this article, we describe the synthesis of light element fully ionized plasma atmosphere models of 
Lloyd (2003) with realistic NS temperature and field configurations derived from NS crustal conductivities 
(Heyl & Hernquist 1998a). The resulting model yields self-consistent pulse properties of the pulsar surface 
radiation and provides substantial improvement over previous efforts to describe pulsed thermal X-rays from 
NSs. We compare the pulsar atmosphere model to archival Chandra data for PSR 0656-1-14, and examine 
the extent to which our model provides an adequate solution to both the average spectral envelope of this 
pulsar and the energy-dependent pulsed fraction. The pulsar model is capable of reproducing the average 
X-ray spectrum without an explicit secondary thermal component, but fails to reconcile the observed pulsed 
fraction. 

We organize our discussion as follows: the pulsar-atmosphere (PA) model is detailed in §2. The ex- 
traction and reduction of the Chandra LETG and ACIS-S observations is described in §3.1, followed by the 
results of our spectral analysis in §4. 



2. Pulsar atmospheres 



The surface distribution of effective temperature on a NS depends on the local magnetic field strength 
and its inclination with respect to the surface normal vector, and is prescribed by the conductive properties 
of the magnetized neutron star crust (e.g. Grcenstein & Hartkc 1983; Hernquist & Applcgatc 1984; Hernquist 
1984; Heyl & Hernquist 1998a). Both the field and temperature distributions are parameterized by their 
values at the magnetic pole. For a body-centered dipolar magnetic field, the surface field strength at magnetic 
colatitude 6^ is 



BiOp) = ^^l + 3cos2ep . 



(1) 



The magnetic field vector at this point is inclined with respect to the local normal (i.e. radial vector) 
by an angle 'ip 



cosV' ■ 



2 cos 9v 



yr+3cos2^ 



(2) 



Thermal conductivity in the magnetized NS crust varies with the magnetic field density and direction, 
resulting in a smooth variation of effective temperature from the magnetic pole to equator; for fields Bp > 
10" G, Heyl & Hernquist (1998a) find: 



1/5 



4 cos^ 0p 
1 + 3 cos^ 6j) 
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(3) 



We may optionally superimpose the radiation from heated polar caps on the pulsar temperature distri- 
bution. The cap dimension is described by the angular radius Oc, and its intensity pattern is modeled with a 
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normal, uniform surface magnetic field Bp and constant effective temperature Tc on the magnetic colatitude 
interval [0,0c]- Radiation from the polar cap will supersede that of the pulsar distribution for Op in this 
range. We simplify integration of the total surface emission in the ordinary dipole field geometry by taking 
advantage of the axial symmetry of the B and T distributions. The radiation field at an arbitrary point on 
the stellar surface can be mapped into the equivalent point on a narrow strip of magnetic longitude extending 
from the (magnetic) pole to equator which we discretize into ~ 20 surface elements having uniquely defined 
local properties B, Tejj from equations (1-3). Each element is modeled using the computational methods 
described in Lloyd (2003) to derive the asymmetric thermal intensity patterns for a hydrogen plasma atmo- 
sphere with these parameters. The model intensities Isik) from all elements are tabulated for interpolation 
of the radiation profile from an arbitrary point on the NS surface. The radiation field is fully symmetric 
about the magnetic equator, including the two antipodal caps. 



2.1. Phase resolved model spectra 

The viewing geometry of the NS is defined by two angles: the angle x between the line of sight (LOS) 
and the pulsar spin axis P, and the relative angular displacement ^ of the magnetic dipole axis B,iip from P. 
Consequently, the angle a between the magnetic dipole axis and the LOS varies with phase angle 7 = 2nt/P: 

Bdip = (sin X cos ^ — cos x sin ^ cos 7, — sin ^ sin 7, cos a) (4) 

where 

cosa = cosxcos^ + sinxsin^cos7. (5) 

Any particular viewing geometry for which x or ^ = is degenerate, and fails to produce modulation of the 
surface emission in the dipole model. Results for choices of x or ^ which exceed 7r/2 are equivalent to those 
for which the angle is subtracted from tt. 

Calculation of the phase-resolved flux proceeds from integration of the intensity vectors coincident with 
a distant observer's LOS from a large number of points on the stellar surface. Two coordinate systems are 
required, and computation of the model spectrum is facilitated by the translation from viewing or integration 
coordinates defined by a fixed LOS to NS surface coordinates defined with respect to the dipole axis. This 
derivation proceeds in Cartesian coordinates and assumes I = (0,0, —1) for the LOS. To deduce the correct 
intensity contribution from an arbitrary point on the star, wc require three surface coordinates: the local 
magnetic colatitude 9p, and the direction of the intensity vector which is parallel to the LOS at distance 
and given by the pair of angles {9b, (pb) in the NS frame. The global geometry is illustrated in Figure 1. The 
remainder of our derivation describes recovery of this angular information from the dipole pulsar model for 
arbitrary x, ^ and 7. 

The precise intensity contribution from each surface element varies with rotational phase for a given 
viewing geometry, but the overall configuration of surface elements remains fixed. In a given phase, each 
surface element is identified by a pair of polar integration coordinates (5, 0) where 5 is the angle between the 
intensity vector from the surface which intersects the detector [0,7r/2] and the local normal vector, and the 
circumpolar angle 4) is on the interval [0, 2-k\. The LOS is the pole of this coordinate system. The extent of 
the visible surface is defined by 5 = 7r/2. Note that, for relativistic stars, photons with trajectories (5 = 7r/2 
in these coordinates emerge from surface latitudes 6 > 7r/2; in such cases the star is self-lensing and has 
larger exposed radiative area than geometric area for a hemisphere: for a non-relativistic star ~~> 5. To 
provide an absolute reference for cj), the prime meridian is defined such that the spin axis P and LOS are 



- 6 - 




Fig. 1. The global viewing geometry is specified by the pair of angles x,^. The dipole axis sweeps the 
dashed curve as the star rotates though phase angle 7. P is fixed and coplanar with the line of sight {—z) 
while Biiip is coplanar with the LOS and spin axis P when 7 = or tt. Any surface element having normal 
h is referenced in two coordinate systems: (a, 6p) for some phase in the pulsar coordinates, and by {6, (j)) in 
the surface integration coordinates. <f) is measured with respect to the projection of h onto the x — y plane 
(shown here as hxy)- 
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coplanar with all normal vectors having (j> = regardless of phase angle; for 7 = 0, the magnetic dipole 
vector is likewise coplanar with the spin and LOS. 

Figure 2 summarizes the geometry at a given surface element. In our integration coordinates, the local 
normal vector at the point {9, (j)) on the stellar surface is 

h = (sin 6 cos (j), sin sin cj), cos 0) (6) 

and has magnetic colatitude given by 

cos dp = I3dip ■ n (7) 

The magnetic field at n is _B = 3 cos Opfi — B^ip (cf. eqn 2). In the atmosphere description of the individual 
surface elements, the local azimuthal reference is = for coplanar n, B. The angle between the LOS and 
n — B plane is 

cos = m-v (8) 

where v = B — hcos^j and m = — (cos^coS(/<, cos sine/), — sinO) arc the projections of B and the LOS onto 
the local tangent plane. is the azimuthal angle the intensity vector forms with respect to the magnetic 
field axis and local normal vectors (Lloyd 2003). Note that the intensity vector and local normal are coplanar 
with the LOS. To complete our definition of the intensity vector, we measure the angle 6b between the desired 
intensity vector at the NS surface and the local normal, given by the general relativistic "ray-tracing" integral 
(Page 1995): 



r-i/2 



0b{x) = X 
Jo 

where x = sin 6, = R/Rs and the Schwarzschild radius Rs = 2GM/(?. 



-1/2 

Av (9) 



The intensity IsiOb, <Pb) from any integration element is step-wise interpolated from the table of surface 
elements. The analytic surface temperature distribution (3) poorly approximates the local flux for dp > 85°; 
this region contributes <C 1% of the total surface emission and integration elements within this interval will 
be considered "dark" making no contribution to the total. The phase resolved spectrum is 

Fe{i) = ^[^)' fj ^[ d(x2) lE{e{x), cj>) (10) 

The prefactors in eqn (10) normalize the total flux to a distant source, where R and d are respectively the 
stellar radius and source distance. Relativistic self-lensing effects are included in the intensity integration 
described above, while the gravitational redshift will be accounted for when fitting the model to the Chandra 
data (§4). The phase averaged spectrum is 

1 Z"^'' 

i^E = ^ I ^iFe{i) (11) 

The average spectrum is typically calculated over 20 — 30 phases of the pulsar rotation 7 = [0, tt] while the 
integration of total flux is performed by interpolating over on the order of 10,000-20,000 surface elements 
per phase. The most useful representation of the phase-resolved spectrum is the energy dependent pulsed 
fraction (PF) 

jpmax ipmin 

PFe = =^ (12) 

^ Tpmax _|_ Tpmin ^ ' 

E E 

In eqn (12) the "bright" (Fg?"^) and "dim" (.F^™) phases are 7 = and 7r/2 respectively which are correct at 
X-ray energies for pulsar models in the absence of a finite polar cap. For cap models with finite temperature 
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Fig. 2. — The surface element defined at the point referenced by 6, (f) in the integration coordinates. The 

intensity / intersects the detector at infinity and the angles 5 and 9 are mutually defined by eqn. (9). 
Individual surface elements are computed using the local B and Teff defined by the global (dipolar) magnetic 
field geometry. The angular intensity patterns are generally anisotropic, and the local azimuthal angle 0b is 
required for interpolation over the emission from the table of computed elements. 
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contrast, the assignment of 7 for bright and dim may be reversed for portions of the X-ray flux (Perna et al. 
2001). This is equivalent to a phase change in the pulse synchronization for different spectral components 
and is distinguished in these models by a negative PF. For a realistic detector the pulsed fraction must be 
evaluated over a finite band: 



Instrumental efficiency and absorption in the ISM are energy dependent and must therefore be included 
in the calculation of the integrated pulsed fraction (Page 1995; Perna et al. 2000). As with gravitational 

redshift of the spectrum, these effects will be incorporated into spectral model from within the XSPEC 
software package. The pulsed fraction predicted in the PA model and others increases with large R, owing 
to the diminished effect of self-lensing. 

For a given polar magnetic field strength, the spectral model is parameterized by the polar temperature 
of the pulsar distribution Tp, the angular size 6c and temperature of the heated cap region, and the 
viewing angles x ^-nd ^. The pulsar atmosphere spectra are written in the FITS table format for use in 
the XSPEC software package (Arnaud 1996). Polar caps may take angular radii up to 13° for temperatures 
log Tc = 5.7 — 6.5, while the underlying pulsar pole temperature takes values log Tp = 5.4 — 6.0. The "bright" 
and "dim" phase spectra are also included in the table to facilitate integration of the pulsed fraction. The 
pulsar viewing angles x, C ni^Y assume values from — 90°. 



PSR 0656+14 was observed with the Chandra HRC-LETG configuration (the combination of the High 
Resolution Camera imaging detector and the Low Energy Transmission Grating - Brinkman et al. 2000) 
on 1999 November 28, with a total exposure time of 38 ks (ObsID 130). We retrieved the primary data 
products of this observation (i.e. standard filtered events files, aspect solution, and bad-time/pixels filters) 
from the public Chandra Data Archive^, and extracted a dispersed spectrum of the source with negative 
ande positive orders co-added, and the background spectrum, using the Ciao (Chandra Interactive Analysis 
of Observations) software (standard shape and size extraction regions were used for both source and back- 
ground). Ancillary Response Functions (ARFs) for negative and positive orders were also built with Ciao, 
and then co-added, while for the redistribution matrix the HRC-LETG combination does not allow us to 
resolve orders, so our final co-added spectrum contains counts from all possible dispersed orders. However, 
contamination from orders higher than the first, is only relevant at wavelengths longer than ~ 80 A (^^ 0.15 
keV), which we do not use in our analysis. For the purpose of spectral analysis, we then built (with Ciao) 
negative and positive 1st order HRC-LETG Ancillary Response Functions (ARFs), and then co-added them. 
For the Redistribution Matrix Function (RMF), instead, we used the "canned" Ist-order HRC-LETG RMF 
distributed with the Chandra CALibration Data Base (CALDB). 

The total number of background subtracted counts in the source dispersed spectrum, between 0.15 and 
1.5 kcV (the energy range used in our spectral analysis) is 9230, corresponding to a signal to noise per 
resolution element of AA = 0.05 A, of SN = 2.5 (6 counts per resolution element). Such a low number of 
counts per resolution element does not allow us to search for narrow absorption or emission features in the 




(13) 



3. 



Data reduction 



3.1. Low energy transmission grating spectrum 



^http://asc.hax vaxd.edu/cda/ 
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spectrum or to apply the statistics to fit models for the continuum emission to the unbinned spectrum. 
We then grouped our spectrum allowing a minimum of 50 counts per new grouped channel, and performed 
spectral fitting by using the fitting package XSPEC on background subtracted data. 

3.2. ACIS-S3 CCD spectrum 

PSR B0656+14 was observed for 25 ks on 15 December 2001, using the ACIS detector in continuous- 
clocking (CC) mode (ObslD 2800). In this mode, one imaging dimension is sacrificed in order to provide fast 
read-out of the CCD. The resulting data consist of one row of events for which each pixel value represents the 
sum of events from the associated CCD column, with a time resolution of ~ 3 ms. The target was centered 
on the S3 chip, which is a backside-illuminated device. The energy resolution is AE/E 0.1 at 1 keV, and 
this varies gradually with row number. 

Source and background spectra were extracted from the data using the dmextract routine in Ciao. 
Source events were extracted from a region extending ~ 1.8 arcsec on either side of the projected position 
of PSR B0656+14, while background events were taken from a ~ 48 arcsec region adjacent to the pulsar. 
Spectra were grouped to contain a minimum of 25 counts in each bin. 

Calibration products for spectra obtained in CC-mode do not currently exist. As a result, we used spec- 
tral response and effective area files created for a Timed- Exposure (TE) mode observation of PSR B0656+14 
carried out on the same day, with the same pointing and detector position configuration, as for the CC-mode 
observation. Spectra from the TE-mode data were not used in these analyses because the count rate for the 
pulsar is sufficiently high (2.07 ± 0.01 ct s~^ in the CC-mode data) that pileup effects significantly distort 
the spectrum. 

3.3. Timing analysis 

For investigation of the pulsed X-ray emission from PSR B0656-I-14 we extracted event times from a 1 

X 1.5 arcsec^ box which included pixels for which the count rate was a factor of 10 or more higher than the 
average for pixels outside the region of the pulsar location. Event times were corrected from readout times 
to arrival times by correcting for time offset effects from dither of the spacecraft and any motions of the 
Science Instrument Module, following the science analysis thread on CC-mode time corrections with Ciao. 
These times were then corrected to the solar system barycenter with the axbary task, using the nominal 
coordinates of the pulsar (RA2000: 06:59:48.122, Dec2ooo: -hl4:14:21.53). 

A search for periodicity around the expected frequency of 2.59805548 Hz, based on the ephemeris of 
Chang & Ho (1999), was carried out using the Z„ test (Buccheri et al. 1983). The pulsations were easily 
detected with / = 2.598054 ± 1.5 x 10~^ Hz. Absolute phases for each event were then calculated based on 
the above ephemeris (see Table 1). Pulse profiles for selected energy ranges are shown in Figure 3, along 
with the measured pulsed fraction for each profile, defined here as p = {M — m)/{M + m) where M is the 
maximum number of counts in the folded lightcurve and m is the minimum. 




Fig. 3. — X-ray lightcurves for PSR 0656+14 accumulated with Chandra ACIS, integrated on representative 
energy intervals. The pulsed fraction for each band is indicated in the upper left corner of the panel. The 
pulse profiles are irregular and are not well described by the high degree of symmetry intrinsic to our pulsar 
atmosphere model. The energy intervals are the same as found in Pavlov et al. (2002). 
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4. Spectral analysis 

Radio measurements of PSR 0656+14 suggest the pulsar has a magnetic field strength B ~ 4.7 x 10^^ G 
in the dipole braking model. Its spectrum has been observed at optical and UV energies, and Ramanamurthy 
et al. (1996) reported the possible detection of 7-ray pulsations at the 385 ms rotation period in the EGRET 
data. Neither atomic (ionization) nor cyclotron features are evident at any wavelength, and a direct measure 
of the surface magnetic field is not possible. For this analysis we assume that the surface polar strength 
does not differ from that estimated from radio measurements, consistent with the absence of either electron 
or proton cyclotron absorption at X-ray energies. Our analysis proceeds from the combined LETG (0.3 — 
0.9 keV) and ACIS-S (0.4 - 6.0 keV) data sets. 

A uniform (single component) blackbody spectrum will not produce modulated flux, and both single 
temperature blackbody (BB) and uniformly magnetized neutron star atmosphere spectra (ATM) also fail to 
reproduce the phase-averaged Chandra spectrum; these elementary models can be discarded immediately. 
The next simplest alternative which does permit modulation is the composite of either uniform thermal 
emission and a power-law (PL) component or a double-blackbody (DBB) configuration with finite tempera- 
ture contrast. The ACIS spectrum is non-thermal above 2.0 keV; the DBB model must also include the PL 
component to obtain a satisfactory fit. We fixed the power-law component to be consistent with the spectral 
fit to optical measurements as described by Kurt ct al. (1998). The ATM-I-PL yields a formally acceptable 
fit to these data, although the resultant power-law amplitude is sufficiently small that temperature variation 
is required to account for the PF below 2.0 keV, even for 100% pulse index in the PL flux at these energies. 
The extent to which the non-thermal emission is pulsed is arbitrary in this model; however, no evidence of 
pulsations above 2 keV was found in RXTE observations of this source (Chang & Ho 1999), and we assume 
the non-thermal flux is unpulsed at all energies in the Chandra data. It is possible that finite pulsation of 
non-thermal emission at energies < 6.0 keV could remain consistent with the RXTE measurement. The 
combined BB+PL is a poor representation of these X-ray data, and we find that although the DBB plus PL 
model provides an acceptable fit to the average spectrum (Ti = 0.85 MK, T2 = 1.5 MK, A2/A1 = 3.4 x 10~^ 
Pavlov et al. 2002), this configuration does not reproduce the observed pulsed fraction at 0.5-1.0 keV for any 
choice of viewing angles. 

We evaluated the combined Chandra data sets with a grid of PA model spectra for a grid of mass 
M — (0.4, 1.0, 1.4)M0, radius i? = 10 — 16 km in steps of 1 km, and temperature log Tp = 5.4 — 6.0 in steps 
of 0.1 dex, and we obtained good spectral fits (x^/ dof < 1.0) with each pair of M, R for a range of viewing 
angle x on the interval 2 — 44°. Here, we have assumed X = C to reduce the parameter space, which maximizes 
the pulsed fraction while demanding the magnetic pole sweep through the line of sight. For each pairing, we 
assumed the non-thermal emission was well represented by a power-law having index 7 = 1.5. The spectral 

Table 1: PSR B0656-I-14 Ephemeris 



Quantity 


Vahie 


Epoch (MJD) 


50701.000002341 


/(Hz) 


2.5981054226747 


/(Hz s-i) 


-3.71150 X 10-13 


/(Hz s-^) 


8.33 X 10-25 
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fits favor low absorbing columns Nh.20 — 1.0—1.2. From this grid, we selected those models which generate 
the observed pulsed fraction from 0.3-0.5 keV, and maximized the source distance (Figure 4). The PA model 
provides a good fit to the observed average spectrum for a range of Tp = 0.48—0.62 MK {T^ = 0.45—0.5 MK) 
yielding a source distance 100 %d< 250 pc. The thermal component PF(0.3 — 0.5 keV) < 0.32 over the grid, 
but also < 0.36 above 1.0 kcV. The deficit of pulsed emission predicted by the PA model may be accounted 
for by finite modulation of the non-thermal radiation at these energies at the level of ~ 70%. 

The distance and pulse index constraints favor sub-solar NS masses and large stellar radii. To proceed, 
we consider in detail the PA model results for a particular choice which satisfies our constraints, M = 0.4M© 
and i? = 16 km. We recovered PF(0.3 — 0.5 keV) within the measured uncertainties in our models for 
Tp ~ 0.5 MK and % = 30 ± 2°. The unfolded spectral model is illustrated in Figure 5. The PA model fiux 
peaks at a lower energy compared to a uniform surface at Tp and has a broader spectral envelope owing 
to the area-weighted average of the surface distribution, allowing for excellent representation of the X-ray 
spectrum. The PA component does not over-estimate the optical flux of the pulsar, which has a non-thermal 
spectrum (Figure 6). 

For X = the ideal piilsc;d fraction predicted by the model exhibits two trends with R and x- Spectral 
fits for fixed R produce larger net pulsed fraction for increasing % owing to larger flux contrast as the bright 
magnetic pole is inclined from the line of sight. This has the incidental effect of requiring slightly larger 
model normalization, favoring lower d as the overall flux is reduced on average. Second, the magnitude of 
sclf-lcnsing is reduced by cither increasing the radius or reducing the stellar mass, thereby increasing the 
pulsed fraction at all energies. The effective emitting area of the star grows with increasing R°° = R{1 + z), 
and spectral fits to the pulsar spectrum accommodate greater source distance. 

The energy dependence of the pulsed fraction (Figure 7) is an especially restrictive constraint on models 
for thermal X-ray production for PSR 0656-1-14. Page (1995) demonstrated how continuous variation of the 
surface temperature naturally produces a PF which increases monotonically with energy. We attempted 
to model the pulsed fraction at all energies by forcing finite polar cap emission with temperature contrast 
Tc> Tp. While formally acceptable fits to the phase-averaged spectrum can be obtained from these model 
combinations, the result in each case demands that the "soft" thermal component have both very low Tp 
and large stellar radius. For example, we recovered 6c — 2° for Tp = 0.4 MK and Tc = 1.0 MK held constant 
during the fit; by evaluating the PF over a grid of (x, £,) found the most compatible values were x = 80°, 
^ = 50°. The PA & cap model is nevertheless unacceptable as a comprehensive representation for the 
pulsar thermal radiation: the radius expected for the relatively cool Tp component is 65 km for d = 200 pc; 
alternatively, the pulsar must have a distance of only 30 pc for the canonical radius R = 10 km. Moreover, 
it is difficult to increase Tp enough to substantially reduce R/d without disturbing the relative component 
amplitudes over 0.5-1.0 keV responsible for the "dip" in pulsation on this interval; there is an effective limit 
Tc < 1.3 MK, above which the heated cap becomes too small to compensate for the observed growth of PF 
with energy. 

We also considered examples of discrete models having uniform temperature and field over all but 
the antipodal caps, which had finite temperature contrast. As for the PA-|-cap models, the temperatures 
and effective areas of the emitting regions were adjusted to reproduce the average spectrum and provide a 
suitable crossover energy. Each of the discrete models requires either an unacceptably large stellar radius 
(i? > 40 km for d = 500 pc), or unreasonably short distance; and none has the desired PF evolution. The 
magnetic field strength of the cap {B = 0, few x 10^^ G) did not affect these results. We find the "hard" 
thermal component cannot have both Teff greater than about 0.8 MK and sufficient amplitude to induce 
the crossover on the required energy interval. 
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160 180 200 220 240 

d (pc) 



Fig. 4. — The grid of PA+PL models, organized by mass and radius. Each curve is the locus of predicted 
PF and d for each choice of M, R and x- 16 km tracks for each mass are labeled, and curves for lower 
radii (in 1 km increments) are displaced to lower d for large inclinations x- For each curve, the viewing 
angle x(= £,) increases from ^ 2° from the bottom of the figure to ^ 35° at the top. WcU-dcfincd spectral 
fits were obtained for 7 = 1.5 and Nh = 1-0 — 1.2 x 10^" cm~^. Distance and PF limits restrict further 
analysis to those models in the rectangular area. The best PA+PL models gives a 15-16 km stellar radius 
for d ~ 190 pc, favors a sub-solar neutron star mass and implies x — 30°. 
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PSR 0656+14 - PA+PL unfolded spectrum 
LETG &c ACIS-S 




^ 0.5 1 2 5 

channel energy (keV) 

Fig. 5. — The unfolded model and spectrum for i? = 16 km enclosing 0.4 Mq at d ~ 190 pc. These stellar 
parameters are consistent with several equations-of-state (Lattimer & Prakash 2001). Non-thermal emission 
dominates above 2 keV and is arbitrarily pulsed in our model. 
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Fig. 6. — The two component optical model for PSR 0656+14. Data points are those compiled in Koptsevich 
et al. (2001). The optical flux is dominated by the non-thermal component, described here by the same 
power-law used to model the hard X-rays in the Chandra data. The thermal component has an approximately 
Rayleigh-Jeans spectrum and does not over-produce optical flux. The broad feature centered at ~ 7 x 10^^ Hz 
is proton cyclotron absorption in the hydrogen plasma. 
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Fig. 7. — The pulsed fraction for PSR 0656+14. The measured net PF integrated from 0.13-6.0 keV is 
13.2%. The convolved model PF for the parameters of Figure 5 (open triangles) fails to reproduce the 
pulsed fraction from 0.5-1.0 kcV. The solid curve illustrates the thermal contribution to the pulsed fraction 
for the same model. Modulation of ~ 70% in the PL component is sufficient to reproduce the measured PF 
above 1 keV (filled triangles), but does not lower the predicted pulsation from 0.5-1.0 keV; the dashed curve 
illustrates the PF{E) for this model. Calibration issues prevent using the spectrum below 0.3 keV, although 
these data are expected to faithfully represent the pulsed fraction at lower energies. 
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5. Discussion 

The amplitude of X-ray pulsations from PSR 0656+14 exclude uniform temperature models for their 
origin, and require either a discrete or continuous distribution of surface temperatures. We have analyzed 
the archival Chandra observations of this middlc-agcd radio pulsar with a variety of dual component thermal 
models and the PA model of §2. We find that the phase- averaged spectrum is well described by the combined 
PA+PL model for a broad range of acceptable M, R with the addition of a power-law component at energies 
> 2 keV. By maximizing the source distance and requiring that the resultant model reproduce the measured 
pulsed fraction at 0.3 — 0.5 keV we find the PA interpretation favors large NS radii i? ~ 15 — 16 km, and 
stellar mass M < 1.0 Mq. These results are compatible with several equations-of-state for the neutron star 
core (Lattimer & Prakash 2001) and also indications from neutron star cooling analyses that the mass of 
PSR 0656 is smaller than the canonical value M = 1.4Mq (e.g. Yakovlev & Haensel 2002; Tsuruta ct al. 
2002). For a particular choice of mass and radius, i? = 16 km, M = 0.4 Mq we recovered the viewing 
angle x = ^ — 30° and a source distance d ~ 190 pc. Serious difficulties remain in reconciling the average 
X-ray spectrum with the measured modulation. In particular, the apparent deficit of pulsed emission at 
0.5 — 1.0 kcV is not readily explained with the PA model without the addition of a heated polar cap; the 
discrete models we investigated fail to generate the correct PF evolution while models based on the pulsar 
type surface can only reproduce the PF for unreasonable (R/d). The model also fails to generate the full 
measure of pulsation observed from 1.0 — 6.0 keV without additional modulation of the non-thermal flux; 
~ 70% in the PL component would be suSicient to close the gap above 1.0 keV. 

The continuum spectrum from an iron envelope or from high metallicity atmospheres bear some resem- 
blance to a blackbody, although such metal rich compositions have pronounced atomic line and ionization 
features which should be evident at X-ray energies (Rajagopal & Romani 1996). These features are not found 
in the data; indeed, the (unmagnetized) Fe models of Gansicke et al. (2002) are a very poor representation 
of the Chandra observations, which favor the light element plasma model we have described. However, the 
presence of strong absorption features in magnetized hydrogen plasma has previously been associated with 
the formation of a soft X-ray crossover energy (Zavlin et al. 1995) above which the pulse profiles are phase 
shifted with respect to those at longer wavelengths, and which also implies a reduced pulse contrast (13). 
For the magnetic field strength 4.7 x 10^^ G, the ground state binding energy of hydrogen is approximately 
0.25 keV (Potekhin 1998), below the usable lower limit of these Chandra data. A magnetic field B ~ 10^^ G 
is required to increase the binding energy (and its putative phase shift) to 0.5 keV; the X-ray data do not 
reveal this feature and the pulsed fraction from 0.5-1.0 keV is not adequately explained by photoionization. 

The lightcurves predicted by the PA model are symmetric about the 7 = phase. The observed X-ray 
lightcurves, however, show irregular pulse profiles which are an indication of the fact that the simple PA 
model is an oversimplification which lacks the necessary variation in surface structure or beaming to generate 
the measured PF at all energies. 

We found a restricted range of viewing angles x — 30° best represented the modulated thermal flux in 
our PA model. This result is reminiscent of estimations for x from very different wavelengths and emitting 
regions for this pulsar (Malov 1990; Rankin 1993); Rankin (1993) also deduced z < 0.15 consistent with the 
present results which favor z ~ 0.04 or smaller. Harding & Muslimov (1998) have estimated the angles x 
and ^ from a curvature radiation model to be smaller than approximately 30°. This is also consistent with 
the estimate x ^ 23.5 — 35° derived by Malov (1990) from three measures of the radio pulse shapes. 

Some models for 7-ray curvature radiation suggest a correlation between the presence of currents 
and the extent to which thermal surface X-rays can be modulated. Models for the influence of gamma-ray 
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radiation mechanisms in the outer magnetospherc provide for the production of X-ray power-law spectra 
and the formation of a blanket of charge through which weakly modulated X-ray surface radiation, and 
more strongly pulsed non-thermal polar cap radiation, would be observed (Wang et al. 1998). Pulsars 
without accelerators for gamma-ray production would be dominated by the surface cooling radiation at 
X-ray energies. PSR 0656+14 is, at best, a weak source of pulsed 7-rays (Ramanamurthy ct al. 1996), 
which does not favor the "hohlraum" picture of Wang et al. (1998); furthermore, the charge blanket scenario 
provides an interpretation for strongly pulsed non-thermal emission, but does not adequately explain the 
breadth of the thermal X-rays. The X-ray model of Cheng & Zhang (1999) describes the pulsar spectrum 
as primarily thermal in nature and produced through heating of the stellar surface by return currents of 
curvature radiation pairs from the polar gap and outer gap. Their model provides a unified mechanism for 
describing both thermal and non- thermal X-ray emissions from PSR 0656-1-14 for a viewing geometry similar 
to our findings. 

An additional problem arises from consideration of the extent to which the atmospheric plasma can be 

heated externally while preserving the thermal character of its net emission. Zanc et al. (2000) have calculated 
spectra from model atmospheres heated externally (e.g. by charged particle bombardment) — the resultant 
X-ray spectra are essentially unaff'ected by heating for modest particle fluence, although amplification of the 
optical spectrum is apparent. We have considered the related problem of illuminated atmospheres, in which 
the surface plasma is subject to an external radiation field of the form Ai/^'^ . The physics of illumination 
differs from that of heating but models for each yield comparable results: heating of the outermost strata of 
the atmosphere with concomitant enhancement of the optical and near UV radiation. The X-ray beaming 
functions (and hence the PF) however are unaffected by the effects of illumination unless the incident flux 
is so great that the atmosphere is completely disrupted. 

It is unlikely that the global magnetic field geometry for a given pulsar is adequately described by the 
centered dipole which is assumed in our pulsar-atmosphere model. Contributions from higher multipole 
moments, displacement of the dipole, stellar oblateness or localized distortions of the field (Page et al. 1995) 
can impart additional pulse components and erode the high degree of symmetry found in the PA results. 
Interaction between thermal radiation and charged clouds in the magnetospherc (Wang et al. 1998), or with 
non-thermal currents may also affect the X-ray lightcurves and deduced PF. 

This work was supported in part by NASA contract NA58-39073 (PS). 
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